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Abstract In this paper we present the solver DuQuad specialized for solving 
general convex quadratic problems arising in many engineering applications. 
When it is difficult to project on the primal feasible set, we use the (augmented) 
Lagrangian relaxation to handle the complicated constraints and then, we ap¬ 
ply dual first order algorithms based on inexact dual gradient information for 
solving the corresponding dual problem. The iteration complexity analysis is 
based on two types of approximate primal solutions: the primal last iterate 
and an average of primal iterates. We provide computational complexity esti¬ 
mates on the primal suboptimality and feasibility violation of the generated 
approximate primal solutions. Then, these algorithms are implemented in the 
programming language C in DuQuad, and optimized for low iteration com¬ 
plexity and low memory footprint. DuQuad has a dynamic Matlab interface 
which make the process of testing, comparing, and analyzing the algorithms 
simple. The algorithms are implemented using only basic arithmetic and logi¬ 
cal operations and are suitable to run on low cost hardware. It is shown that if 
an approximate solution is sufficient for a given application, there exists prob¬ 
lems where some of the implemented algorithms obtain the solution faster than 
state-of-the-art commercial solvers. 
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1 Introduction 

Nowadays, many engineering applications can be posed as convex quadratic 
problems (QP). Several important applications that can be modeled in this 
framework such us model predictive control for a dynamical linear system 
an d its dual called moving horizon estimation jTDj, DC op¬ 
timal power flow problem for a power system m, linear inverse problems 
arising in many branches of science nnm or network utility maximization 
problems |30j have attracted great attention lately. Since the computational 
power has increased by many orders in the last decade, highly efficient and 
reliable numerical optimization algorithms have been developed for solving the 
optimization problems arising from these applications in very short time. For 
example, these hardware and numerical recent advances made it possible to 
solve linear predictive control problems of nontrivial sizes within the range 
of microseconds and even on hardware platforms with limited computational 
power and memory m- 

The theoretical foundation of quadratic programming dates back to the work 
by Frank & Wolfe [BJ. After the publication of the paper [B] many numerical 
algorithms have been developed in the literature that exploit efficiently the 
structure arising in this class of problems. Basically, we can identify three 
popular classes of algorithms to solve quadratic programs: active set methods, 
interior point methods and (dual) first order methods. 

Active set methods are based on the observation that quadratic problems with 
equality constraints are equivalent to solving a linear system. Thus, the itera¬ 
tions in these methods are based on solving a linear system and updating the 
active set (the term active set refers to the subset of constraints that are sat¬ 
isfied as equalities by the current estimate of the solution). Active set general 
purpose solvers are adequate for small-to-medium scale quadratic problems, 
since the numerical complexity per iteration is cubic in the dimension of the 
problem. Matlab’s quadprog function implements a primal active set method. 
Dual active set methods are available in the codes mm- 

Interior point methods remove the inequality constraints from the problem 
formulation using a barrier term in the objective function for penalizing the 
constraint violations. Usually a logarithmic barrier terms is used and the re¬ 
sulting equality constrained nonlinear convex problem is solved by the Newton 
method. Since the iteration complexity grows also cubically with the dimen¬ 
sion, interior-point solvers are also the standard for small-to-medium scale 
QPs. However, structure exploiting interior point solvers have been also de¬ 
veloped for particular large-scale applications: e.g. several solvers exploit the 
sparse structure of the quadratic problem arising in predictive control (CVX- 
GEN [IT], FORCES [5]). A parallel interior point code that exploits special 
structures in the Hessian of large-scale structured quadratic programs have 
been developed in 0. 

First order methods use only gradient information at each iterate by com¬ 
puting a step towards the solution of the unconstrained problem and then 
projecting this step onto the feasible set. Augmented Lagrangian algorithms 
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for solving general nonconvex problems are presented in the software pack¬ 
age Lancelot [5]. For convex QPs with simple constraints we can use primal 
first order methods for solving the quadratic program as in [28]. In this case 
the main computational effort per iteration consists of a matrix-vector prod¬ 
uct. When the projection on the primal feasible set is hard to compute, an 
alternative to primal first order methods is to use the Lagrangian relaxation 
to handle the complicated constraints and then to apply dual first order algo¬ 
rithms for solving the dual. The computational complexity certification of first 
order methods for solving the (augmented) Lagrangian dual of general convex 
problems is studied e.g. in [immm and of quadratic problems is stud¬ 
ied in izrag. In these methods the main computational effort consists of 
solving at each iteration a Lagrangian QP problem with simple constraints for 
a given multiplier, which allows us to determine the value of the dual gradient 
for that multiplier, and then update the dual variables using matrix-vector 
products. For example, the toolbox FiOrdOs [25] auto-generates code for pri¬ 
mal or dual fast gradient methods as proposed in [23] ■ The algorithm in [23] 
dualizes only the inequality constraints of the QP and assumes available a 
solver for linear systems that is able to solve the Lagrangian inner problem. 
However, both implementations [23j[28] do not consider the important aspect 
that the Lagrangian inner problem cannot be solved exactly in practice. The 
effect of inexact computations in dual gradient values on the convergence of 
dual first order methods has been analyzed in detail in [T51HB1H5] . Moreover, 
most of these papers generate approximate primal solutions through averag¬ 
ing [15][l6l[20] . On the other hand, in practice usually the last primal iterate 
is employed, since in practice these methods converge faster in the primal last 
iterate than in a primal average sequence. These issues motivate our work here. 
Contributions. In this paper we analyze the computational complexity of sev¬ 
eral (augmented) dual first order methods implemented in DuQuad for solving 
convex quadratic problems. Contrary to most of the results from the litera¬ 
ture 1213122 ], our approach allows us to use inexact dual gradient information 
(i.e. it allows to solve the (augmented) Lagrangian inner problem approxi¬ 
mately) and therefore is able to tackle more general quadratic convex problems 
and to solve practical applications. Another important feature of our approach 
is that we provide also complexity results for the primal latest iterate, while 
in much of the previous literature convergence rates in an average of primal 
iterates are given |151ll6ll20ll23j . We derive in a unified framework the compu¬ 
tational complexity of the dual and augmented dual (fast) gradient methods 
in terms of primal suboptinrality and feasibility violation using inexact dual 
gradients and two types of approximate primal solutions: the last primal iter¬ 
ate and an average of primal iterates. From our knowledge this paper is the 
first where both approaches, dual and augmented dual first order methods, 
are analyzed uniformly. These algorithms are also implemented in the effi¬ 
cient programming language C in DuQuad, and optimized for low iteration 
complexity and low memory footprint. The toolbox has a dynamic Matlab 
interface which make the process of testing, comparing, and analyzing the al¬ 
gorithms simple. The algorithms are implemented using only basic arithmetic 
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and logical operations and thus are suitable to run on low cost hardware. The 
main computational bottleneck in the methods implemented in DuQuad is the 
matrix-vector product. Therefore, this toolbox can be used for solving either 
QPs on hardware with limited resources or sparse QPs with large dimension. 
Contents. The paper is organized as follows. In section [2] we describe the 
optimization problem that we solve in DuQuad. In Section [3] we describe the 
the main theoretical aspects that DuQuad is based on, while in Section [6] we 
present some numerical results obtained with DuQuad. 

Notation. For x, y G R” denote the scalar product by (x, y) = x T y and the 
Euclidean norm by ||x|| = y/x T x. Further, [it]x denotes the projection of u onto 
convex set A' and distx(w) = ||u— [u]x|| its distance. For a matrix G € R mX11 
we use the notation ||G|| for the spectral norm. 


2 Problem formulation 

In DuQuad we consider a general convex quadratic problem (QP) in the form: 

F* = minF(u) ( = \u T Qu + q T u ) (1) 

ubu V 2 ) 

s.t: Gu T g tE /C, 

where F : R" —> R is a convex quadratic function with the Hessian Q H, 
G £ R pxra , U CM" is a simple compact convex set, i.e. a box U = [lb ub], and 
K, is either the cone R p or the cone {0}. Note that our formulation allows to 
incorporate in the QP either linear inequality constraints K. = R p (arising e.g. 
in sparse formulation of predictive control and network utility maximization) 
or linear equality constraints K. = {0} (arising e.g. in condensed formulation 
of predictive control and DC optimal power flow). In fact the user can define 
linear constraints of the form: lb < Gu + g < ub and depending on the values 
for lb and ub we have linear inequalities or equalities. Throughout the paper 
we assume that there exists a finite optimal Lagrange multiplier A* for the QP 
© and it is difficult to project on the feasible set of problem ©: 

Xqp = {u € U : Gu + g € /C}. 

Therefore, solving the primal problem © approximately with primal first 
order methods is numerically difficult and thus we usually use (augmented) 
dual first order methods for finding an approximate solution for ©. By moving 
the complicating constraints Gu + g € K, into the cost via Lagrange multipliers 
we define the (augmented) dual function: 

dp(\) = min£ p (u, A), (2) 

nGC/ 

where C p (u,X) denotes the (augmented) Lagrangian w.r.t. the complicating 
constraints Gu + g € /C, i.e.: 

C p (u, A) = min F(u) + (A, Gu + g — s) + ^|| Gu + g — s|| 2 (3) 

sE/C 2 
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where the regularization parameter p > 0. We denote s(u, A) = arg min(A, Gu-\- 

s^K. 

g — s) + |||Gu + g — s|| 2 and observe that: 


s(u, A) 


Gu + g + A A 

0 


if p > 0 
if p = 0. 


Using this observation in the formulation d3|), we obtain: 

C p (u,X)= /^( u ) +f dist 'c(G'n + ff+iA) -^||A|| 2 , if p > 0 ^ 

|f(u) + (A, Gu + g), ifp = 0. 

In order to tackle general convex quadratic programs, in DuQuad we consider 
the following two options: 

Case 1: if Q >~ 0, i.e. Q has the smallest eigenvalue A m i n (Q) > 0, then we 
consider p = 0 and recover the ordinary Lagrangian function. 

Case 2: if Q y 0, i.e. Q has the smallest eigenvalue A m i n (Q) = 0, then we 
consider p > 0 and recover the augmented Lagrangian function. 

Our formulation of the (augmented) Lagrangian (U) and the previous two 
cases allow us to thereat in a unified framework both approaches, dual and 
augmented dual first order methods, for general convex QPs. We denote by 
u( A) the optimal solution of the inner problem with simple constraints u £ U: 

u(A) = arg min £ p (it, A). (5) 

u£U 


Note that for both cases described above the (augmented) dual function is dif¬ 
ferentiable everywhere. Moreover, the gradient of the (augmented) dual func¬ 
tion d p ( A) is Ld-Lipschitz continuous and given by [T51[TB1I^I1E51 : 

Vdp(A) = Gu( A) + 5 - s(u( A), A) and L d = ( 6 ) 

for all A £ R p . Since the dual function has Lipschitz continuous gradient, 
we can derive bounds on d p in terms of a linear and a quadratic model (the 
so-called descent lemma) [2Tj: 

0< [d p (g,) + (X7d p (g),\~g)}-d p (\) < ^f\\p - A\\ 2 VA^eF. (7) 

Descent lemma is essential in proving convergence rate for first order methods 
m- Since we assume the existence of a finite optimal Lagrange multiplier for 
©. strong duality holds and thus the outer problem is smooth and satisfies: 


where 


F* = max dpi A), 
Ae/Cd 


IC d 


R+, if Q F 0 and K = R p _ 
K. p , otherwise. 


(8) 






6 


Ion Necoara, Andrei Patrascu 


Note that, in general, the smooth (augmented) dual problem (JS]) is not a QP, 
but has simple constraints. We denote a primal optimal solution by u* and 
a dual optimal solution by A*. We introduce A* C /C d as the set of optimal 
solutions of the smooth dual problem ([5]) and define for some A 0 £ R p the 
following finite quantity: 

Tl d = min ||A*-A°||. (9) 

A*eA* 

In the next section we present a general first order algorithm for convex opti¬ 
mization with simple constraints that is used frequently in our toolbox. 


2.1 First order methods 


In this section we present a framework for first order methods generating an 
approximate solution for a smooth convex problem with simple constraints in 
the form: 


= min 
x£X 


( 10 ) 


where cf> : R™ — > R is a convex function and X is a simple convex set (i.e. the 
projection on this set is easy). Additionally, we assume that <j> has Lipschitz 
continuous gradient with constant L^ > 0 and is strongly convex with constant 
<T 0 > 0. This general framework covers important particular algorithms [2ll21| : 
e.g. gradient algorithm, fast gradient algorithm for smooth problems, or fast 
gradient algorithm for problems with smooth and strongly convex objective 
function. Thus, we will analyze the iteration complexity of the following general 
first order method that updates two sequences ( x k ,y k )k>o as follows: 


Algorithm FOM (0, X) 

Given x° = y 1 £ X , for k > 1 compute: 


1 . x k = 


y k - T~XHv k ) 


J A 


2. y k+1 = x k + /3k{x k — X k 1 ), 


where /3k is the parameter of the method and we choose it in an appropriate 
way depending on the properties of function (f>. More precisely, /3fc can be 
updated as follows: 

GM: in the Gradient Method Bk = B „ k ~ 1 . where 9^ = 1 for all k. This is 
equivalent with f3k = 0 for all k. In this case y k+1 = x k and thus we have the 
classical gradient update: x k+1 = [ x k — j^X(j){x k )\x- 

FGM: in the Fast Gradient Method for smooth convex problems /?& = e 8 k * , 

where 9k +i = 1 + ^ 2 +4 ~ and 9\ = 1. In this case we get a particular version 
of Nesterov’s accelerated scheme m that updates two sequences (x k ,y k ) and 
has been analyzed in detail in [ 2 ]. 
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fgm ct : in fast gradient algorithm for smooth convex problems with strongly 

convex objective function, with constant > 0, we choose fik = 

for all k. In this case we get a particular version of Nesterov’s accelerated 
scheme m that also updates two sequences ( x k ,y k ). 

The convergence rate of Algorithm FOM(^,X) in terms of function values is 
given in the next lemma: 




Lemma 1 W\2tf For smooth convex problem ©3 assume that the objective 
function cj) is strongly convex with constant > 0 and has Lipschitz continu¬ 
ous gradient with constant > 0. Then, the sequences (a? fc , j/ fc ) fc>0 generated 
by Algorithm FOM^,X ) satisfy: 


(f{x k ) 


cf>* < min 



2 L+Kl \ 

(k+l)p(M J 


( 11 ) 


where 7= min^ ||i° — x*A, with Xf the optimal set of m, and p(/3k ) is 
defined as follows: 


Pifik) 


1 ifPk=0 

2 otherwise. 


( 12 ) 


Thus, Algorithm FOM has linear convergence provided that > 0. Other¬ 
wise, it has sublinear convergence. 


3 Inexact (augmented) dual first order methods 

In this section we describe an inexact dual (augmented) first order framework 
implemented in DuQuad, a solver able to find an approximate solution for the 
quadratic program ©• For a given accuracy e > 0, u e £ U is called an e-primal 
solution for problem © if the following inequalities hold: 

dist ic(Gu e + g) < 0(e ) and | F(u e ) — F*| < 0(e). 

The main function in DuQuad is the one implementing the general Algorithm 
FOM. Note that if the feasible set Xqp of © is simple, then we can call 
directly FOM (F,Xqp) in order to obtain an approximate solution for ©. 
However, in general the projection on Xqp is as difficult as solving the origi¬ 
nal problem. In this case we resort to the (augmented) dual formulation © for 
finding an e-primal solution for the original QP ©. The main idea in DuQuad 
is based on the following observation: from © © we observe that for com¬ 
puting the gradient value of the dual function in some multiplier A, we need to 
solve exactly the inner problem ©; despite the fact that, in some cases, the 
(augmented) Lagrangian C p (u , A) is quadratic and the feasible set U is simple 
in ©, this inner problem generally cannot be solved exactly. Therefore, the 
main iteration in DuQuad consists of two steps: 
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Step 1 : for a given inner accuracy ei n and a multiplier g £ R p solve approx¬ 
imately the inner problem © with accuracy e; n to obtain an approximate 
solution u(g) instead of the exact solution u(g), i.e.: 

0 < Cp(u(g), g) - d p (g) < e in . (13) 

In DuQuad, we obtain an approximate solution u{g) using the Algorithm 
FOM(£ p (’,/z),f/). From, Lemma © we can estimate tightly the number of 
iterations that we need to perform in order to get an ei n -solution u(g) for 
© : the Lipschitz constant is Lc = A max (<3) + p||G|| 2 , the strong convexity 
constant is ac = A m i„(<3 + pG T G) (provided that e.g. 1C = {0}) and 7 Zc < Du 
(the diameter of the box set U). Then, the number of iterations that we need 
to perform for computing u{g) satisfying (1131) can be obtained from dill) . 

Step 2: Once an £i n -solution u(g) for ([5]) was found, we update at the outer 
stage the Lagrange multipliers using again Algorithm FOM(d p , Kd)- Note that 
for updating the Lagrange multipliers we use instead of the true value of the 
dual gradient Vd p (/z) = Gu{g) + g — s(u(g), g), an approximate value given 
by: Wdp(g) = Gu{g) + g- s{u(g), g). 

In SmSCESlIIHl it has been proved separately, for dual and augmented dual 
first order methods, that using an appropriate value for ej n (depending on the 
desired accuracy e that we want to solve the QP 0) we can still preserve the 
convergence rates of Algorithm FOM(d p , JCd ) given in Lemma© although we 
use inexact dual gradients. In the sequel, we derive in a unified framework 
the computational complexity of the dual and augmented dual (fast) gradient 
methods. From our knowledge, this is the first time when both approaches, 
dual and augmented dual first order methods, are analyzed uniformly. First, 
we show that by introducing inexact values for the dual function and for its 
gradient given by the following expressions: 

d p {g) = C p (u(g),g) and X7d p (g) = Gu(g) + g - s(u(g), g), (14) 

then we have a similar descent relation as in 0 given in the following lemma: 

Lemma 2 Given e ln > 0 such that (in holds, then based on the definitions 
m we get the following inequalities: 

0 < [d P (g)+ (Vd p (/z), A - g)} - d p (X) <L d \\g- A|| 2 +2e, n VA, geW. (15) 
Proof From the definition of dp, cn and (fTH) it can be derived: 
dp{ A) = min F(u) + (A, Gu + g - s) + ^\\Gu + g - s|| 2 

net/, se/c z 

< F(u(g)) + (A ,Gu(g)+g- s(u(g),g)} + | \\Gu(g) + g - s(u(g), g)\\ 2 

= F p {u{g),g) + ( Gu{g ) + g - s{u{g),g), A - g) 

= d p (g) + (Vdp(g), A - g), 

which proves the first inequality. In order to prove the second inequality, let 
u £ U be a fixed primal point such that C p (u, g) > d p {g). Then, we note that 
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the nonnegative function h(g) = C p (u,/j,) — d p (fd) > 0 has Lipschitz gradient 
with constant Ld and thus we have EH: 

2 ^ II (Gw + g — s(u, /x)) — Vd p (/z)|| 2 = ^H| V/i(/x)|| 2 
< h(u) — min h(is) < C 0 (u, u) — d 0 (u). 

Taking now u = u{fi) and using (fl3l) . then we obtain: 

IIVd p (/x) - Vd p (/x)|| < \/2LdCi n V/x € K p . (16) 

Furthermore, combining (USD with 0 and (DSD we have: 

rfp(A) > d p (/x) + <Vdp(/x), A - /x) - ^||A - /x|| 2 - e in 

> d p (n) + (Vdp(p), A - /x) - ^||A- /x|| 2 + (Vd p (/x) - Vd p (fj,), A - /x) - e in 

> dp(id) + (X7d p (n), A - /x> - ^||A- /x|| 2 - ||W p (/x)-WpGx)||||A - /x|| - e in 

GeJ r, : _ 

> d P (/x) + (Vdp(n), A - /x) - —1|A - /x|| 2 - 1 /2L d e in ||A - /x|| - e in . 

Using the relation \/a& < (a + b )/2 we have: 

dp(A) > dp(Ax) + (Vd p (/x), A - /x) - T d ||A - /x|| 2 - 2e in , 
which shows the second inequality of our lemma. □ 


This lemma will play a major role in proving rate of convergence for the 
methods presented in this paper. Note that in (1151) ej n enters linearly, while 
in nun e; n enters quadratically in the context of augmented Lagrangian and 
thus in the sequel we will get better convergence estimates than those in the 
previous papers. In conclusion, for solving the dual problem 0 in DuQuad 
we use the following inexact (augmented) dual first order algorithm: 


Algorithm DFOM (d p ,JCd) 

Given A 0 = g} € ICd, for k > 1 compute: 

1. u k satisfying (fTdl) for g = g k , i.e. u k = u(g k ) 


2. X k = 


^ + 2 T-^dM) 


J/c d 


3. /+ 1 = X k + Pk{X k - A fe_1 ). 


Recall that u k = u(g k ) satisfying the inner criterion (fT3l) and Vd p (/x fe ) = 
Gu k + g — s(u k , g k ). Moreover, /3k is chosen as follows: 

DGM: in (augmented) Dual Gradient Method /3k = where 9k = 1 for 
all fc, or equivalently /3k = 0 for all k , i.e. the ordinary gradient algorithm. 
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- DFGM: in (augmented) Dual Fast Gradient Method /3k = where 

i -i- /TTW 

Ok- i-i = —Ai_- - an( i = i- e - a variant of Nesterov’s accelerated 

scheme. 

Therefore, in DuQuad we can solve the smooth (augmented) dual problem © 
either with dual gradient method DGM (/3k = 0) or with dual fast gradi¬ 
ent method DFGM (/3k is updated based on Ok). Recall that for computing 
u k in DuQuad we use Algorithm FOM(£ p (-, pf), U) (see the discussion of 
Step 1). When applied to inner subproblem ©, Algorithm FOM (C p (-, pf),U) 
will converge linearly provided that <jc > 0. Moreover, when applying Algo¬ 
rithm FOM(£ p (-, p k ), U ) we use warm start: i.e. we start our iteration from 
previous computed Combining the inexact descent relation (fT5l) with 

Lemma [1] we obtain the following convergence rate for the general Algorithm 
DFOM(d p ,/Cd) in terms of dual function values of ©: 

Theorem 1 (Z[ H<5UT£| / For the smooth (augmented) dual problem © the 

dual sequences (X k , P- k ) k>0 generated by Algorithm DFOM(d p , Kd) satisfy the 
following convergence estimate on dual suboptimality: 

+ (ij) 

where recall IZd = min ||A* — A°|| and p(/3k) is defined as in m- □ 


Note that in 21 Theorem 2], the convergence rate of DGM scheme is provided 

A k 

in the average dual iterate X k = -jT^y Y) ^ an d not in the last dual iterate 

3=0 

X k . However, for a uniform treatment in Theorem Q] we redefine the dual final 
point (the dual last iterate X k when some stopping criterion is satisfied) as 
follows: X k = 


A fc +2£lVdp(A fe ) 


IC d 


3.1 How to choose inner accuracy ei n in DuQuad 


We now show how to choose the inner accuracy ej n in DuQuad. From Theorem 
[Qwe conclude that in order to get e-dual suboptimality, i.e. F* — d p (X k ) < e, 
the inner accuracy e\ n and the number of outer iteration fc ou t (i.e. number of 
updates of Lagrange multipliers) have to be chosen as follows: 


£in 



if DGM 
if DFGM, 


^^5. if DGM 

e 

if DFGM. 


(18) 


Indeed, by enforcing each term of the right hand side of m to be smaller 
than | we obtain first the bound on the number of the outer iterations £) out . By 
replacing this bound into the expression of e; n , we also obtain how to choose 
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ein, i.e the estimates (fTHl) . We conclude that the inner QP (0 has to be solved 
with higher accuracy in dual fast gradient algorithm DFGM than in dual 
gradient algorithm DGM. This shows that dual gradient algorithm DGM 
is robust to inexact (augmented) dual first order information, while dual fast 
gradient algorithm DFGM is sensitive to inexact computations (see also Fig. 
[[]). In DuQuad the user can choose either Algorithm DFGM or Algorithm 
DGM for solving the (augmented) dual problem 0 and he can also choose 
the inner accuracy for solving the inner problem (in the toolbox the default 
values for ei n are taken of the same order as in (fTHll ). 



Fig. 1: Behavior of Algorithms DGM (left), DFGM(right) in terms of primal 
suboptimality w.r.t. inner accuracy for a strongly convex QP with e = 0.01. 


4 How to recover an e-primal solution in DuQuad 


It is natural to investigate how to recover an e-primal solution for the original 
QP 0. Since dual suboptimality is given in the last dual iterate X k , it is 
natural to consider as an approximate primal solution the last primal iterate 
generated by Algorithm DFOM in X k , i.e.: 

u k = u( X k ). (19) 

Note that the last primal iterate u k = u(X k ) coincides with u k = u(n k ) for 
Algorithm DGM. However, for Algorithm DFGM these two sequences are 
different, i.e. u k ^ u k . We will show below that the last primal iterate u k is 
an -^e-primal solution for the original QP 0, provided that F* — d p {X k ) < e. 
We can also construct an approximate primal solution based on an average of 
all previous primal iterates generated by Algorithm DFOM , i.e.: 


k 

£ 

i=i 


Sk 


Sk=Y.°:r 

i=i 


(20) 
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Recall that Oj = 1 in Algorithm DGM and 9j is updated according to the 

l_j_ /X | 

rule 6 j+1 = v 2 1 and = 1 in Algorithm DFGM. In the sequel, we 
prove that the average of primal iterates sequence u k is an e-primal solution 
for the original QP |T]), provided that F* — d p (X k ) < e. 

Before proving primal rate of convergence for Algorithm DFOM we derive a 
bound on ||A fe+1 — A*||, with X k generated by algorithm DFOM, bound that 
will be used in the proofs of our convergence results. In the case of DGM, 
using its particular iteration, for any A € we have: 

||A fc+1 - A|| 2 = ||A fc - A|| 2 + 2{X k+1 - X k , X k - A) + ||A fe+1 - A fc || 2 

= ||A fc - A|| 2 + 2(A fc+1 - A fc , A fc+1 - A) - ||A fe+1 - A fe || 2 

< ||A fc - A|| 2 + d-(Vd p ( X k ), X k+1 - A) - ||A fc+1 - A fe || 2 

Ld 

= \\X k -X\\ 2 + ^-(Wd p (X k ),X k -X) (21) 

+ ((Vd p (A fc ), A fc+1 - A fc ) - L d ||A fc+1 - A fc || 2 ) 

<||A fe -A|| 2 + ^(d p (A fc+1 )-d p (A)) + ^. 


Taking now A = A* and using an inductive argument, we get: 


II A fc 


A* || < Rd, + 



( 22 ) 


On the other hand, for the scheme DFGM, we introduce the notation l k = 
A fe_1 + Ok(X k — A fc_1 ) and present an auxiliary result: 


Lemma 3 

with 9k+i = 


Let (X k , p k ) be generated by Algorithm DFOM(d pi ICd) 
1+ ^ 2 +46> ~ > then for any Lagrange multiplier X £ we have: 


k 


Ol(d P ( A) - d P (A fc )) + > : 9 j A{\,p*)+L d \\l k - A|| 2 <L d ||A° - A|| 2 + 2^ 0 2 e 

i=o 


3=0 


for all k > 0, where A( A, p) = d p (p) + (Vd p (//), A — p) — d p ( A). 

Using this result and a similar reasoning as in [19] we obtain the same relation 
(E!H) for the scheme DFGM. Moreover, for simplicity, in the sequel we also 
assume A 0 = 0. In the next two sections we derive rate of convergence results 
of Algorithm DFOM in both primal sequences, the primal last iterate (THU) 
and an average of primal iterates m- 
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4.1 The \Je convergence in the last primal iterate u k 

In this section we present rate of convergence results for the Algorithm DFOM, 
in terms of primal suboptimality and infeasibility for the last primal iterate 
u k defined in m, provided that the relations (flSl) hold. 

Theorem 2 Let e > 0 be some desired accuracy and u k = u(X k ) be the primal 
last iterate sequence generated by Algorithm DFOM(d p ,lCd) using the inner 
accuracy from m- Then, after k ou t number of outer iterations given in GSD. 
u kout is y/e-primal optimal for the original QP (flj). 

Proof Using the Lipschitz property of the gradient of d p (-), it is known that 
the following inequality holds m- 

dpi A) < d p (p) + (Vd p (p), A — p) — -L|| Vd p (A) - Vd»|| 2 VA, p £ 

Taking p = A* and using the optimality condition (Vd p ( X*),p — A*} < 0 for 
all p £ ICd, we further have: 

||Vd p (A) - Vd p (A*)|| < y/2L d (F* - dp( A)) VA £ K d . (23) 

Considering X = X k and observing that s(u(X k ), X k ) + Gu* + g £ K, we obtain 
a link between the primal feasibility and dual suboptimality: 

d K (Gu k +g)< \\Gu(X k ) + g- s{u(X k ), X k ) - Gu* - g\\ 

= ||Vd p (A fc )-Vd p (A*)|| 

< ||Vd p (A fe ) - Vd p (A fe )|| + ||Vd p (A fc ) - Vd p (A*)|| 

HU+dToJ _ /- 

< V2L d e in + ^/2L d (F*-d p (X k )). (24) 

Provided that F* — d p {X kout ) < e and using e; n as in (fTHl> . we obtain: 

d K {Gu k r' + 9 )< max j ^^e 3/4 | + (2 L d efl\ (25) 

Secondly, we find a link between the primal and dual suboptimality. Indeed, 
we have for all A £ K, d '- 

F* = min F(u) + (A*, Gu + g — s) 

u€lU,s€lKL 

<F(u(X k )) + (X*,Gu(X k ) + g- [Gu(A fc ) + g] K ) . (26) 

Further, using the Cauchy-Schwartz inequality, we derive: 

F{u k o”‘) - F* > — 11 A* 11 distjc (Gm( A feout ) + g) 

> _^ max |&£lh,_l4_ e 3,4j _ Rd(2 ije) l/3. (2J) 
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On the other hand, from the concavity of d p {-) we obtain: 

F{u{X k )) — F* < d p (X k ) — F* — (Vd p (A fc ), A fe ) 

< d p (X k ) — F* — (Vd p (A*), A fe ) + (Vd p (A*) - Vd p (A fe ), A fc ) + e in 

< d p (X k ) — F* — (Vd p (A*), A fc - A*} + ||Vd p (A*) - Vrf p (A fe )||||A fc || + e in 

< l|A* : ||||Vdp(A fc ) — Vd p (A*)|| + e in 

{24j _ i - 

< \\X k \\^2L^~+\\X k \\^2L d (F*-d p {X k )) + e in . (28) 

Taking k = k ou t and e; n as in (fT8l) , based on (1771) and on the implicit assump¬ 
tion that fc out > 1, we observe that ||A feout || < ||A fcout — A*|| + ||A*|| < 4 R d for 
both schemes DGM and DFGM. Therefore, (12811 implies: 

F(u kont ) — F* ? AR d max | ^=—■’ ^y /2 ^ + 4 ^( 2 ^e) 1/2 + e in . 

As a conclusion, from (1271) and the previous inequality, we get the bound: 

I F(u k ^) - F *| < 4i? d max (3 ^ )1/2 e 3/4 | + 4i? d (2L d e) 1 /2 + Cin , 

which implies |F(zi* out ) — F*| < G(-^/e). Using this fact and the feasibility 
bound (1751) . which also implies distjc(GWg° ut +g) < O(^e), we finally conclude 
that the last primal iterate u kout is -^/e-primal optimal. □ 

We can also prove linear convergence for algorithm DFOM provided that 
Amin(Q) > 0 (i.e. the objective function is smooth and strongly convex) and 
U = R” (i.e. the inner problem is unconstrained). In this case we can show 
that the dual problem satisfies an error bound property mm- Under these 
settings DFOM is converging linearly (see nzmiEa for more details). 


4.2 The e convergence in the average of primal iterates u k 

Further, we analyze the convergence of the algorithmic framework DFOM 
in the average of primal iterates ii k defined in (1201) . Since we consider differ¬ 
ent primal average iterates for the schemes DGM and DFGM, we analyze 
separately the convergence of these methods in u k . 

Theorem 3 Let e > 0 be some desired accuracy and u k be the primal av¬ 
erage iterate given in CUD, generated by algorithm DGM, i.e. Algorithm 
DFOM(d p ,K, d ) with 9^ = 1 for all k > 0, using the inner accuracy from (1181) . 
Then, after k out number of outer iterations given in (US. u kout is e-primal 
optimal for the original QP ©. 
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Proof First, we derive sublinear estimates for primal infeasibility for the av- 

k 

erage primal sequence u k (recall that in this case u k = prey 77,7 ) • Given the 

j -o 

definition of A - 7 " 1-1 in Algorithm DFOM(d p , K, d ) with 9j = 1, we get: 


A 7+1 = 


A ' + ^kf dp{XJ) 


Vi > 0 . 


-I/Cd 


Subtracting A - 7 from both sides, adding up the above inequality for j = 0 to 
j = fc, we obtain: 


1 


3 = 0 L 


A j + Vdp(A J ) 


- A J 


J/c d 


k + 1 


||A fc+1 — A 0 1|. (29) 


If we denote z - 7 = A - 7 + 57 —Vd^A- 7 ) — A - 7 + 2 j-Vd p (A J ) 


, then we observe 


that z 3 G /C. Thus, we have X] ^ € /C. Using the definition of Vd p (A- 7 ), 


we obtain: 


j=o 


distAc(Gu^ + g) < 


Ai A 

X £(g»> + 9 ) - fti S ( 2i ' 5 '+*(«’'■ ^)) 

3=0 


k + 1 

i-^(W p (A J )- 2 V) 


i=o 

k 


k + 1 


3=0 


1231 Sid IIAfc+1 0| 


k + 1 


||A fc+1 — A° ||. 


Using ||A fc — A 0 1| < ||A fc — A*|| + Rd and the bound (l22l) for the values e; n and 
k = fc out from (TT51) in the previous inequality, we get: 


dist/c (Gu kout + g) < 


4L d i?d 


TdCin e 

^out Rd 


(30) 


It remains to estimate the primal suboptimality. First, to bound below F(il kout )— 
F* we proceed as follows: 


F* = min F(u) + (A*, Gu + g — s) 

u£U,s£JC 

< Hu 1 :) + {A*,Gu k + g - [Gu k + g] K ) 
<F(u k e ) + \\A*\\\\Gu k e +g-\Gu k e +g\ lc \\ 
= F(u k ) + R d distjc (Gu k + g). 


Combining the last inequality with (13011 . we obtain: 

- e < F(u kont ) - F*. 


(31) 
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Secondly, we observe the following facts: for any u € U, d p ( A) < F* and the 
following identity holds: 

d P ( A) - (Vdp(A), A) = F(u( A)) + |||Vd p (A)|| 2 > F{u( A)). (32) 

Based on previous discussion, m and (1521) . we derive that 
||A fe+1 - A || 2 

? l|A fe - A|| 2 + A- {dpiX k+1 ) - d(X k ) + <Vd p (A fc ), A fc - A) + e in ) 

Ad 

lf32l 1 / n _ x 

< HA fe - A|| 2 + — (F* - F(u k ) - |||W p (A fe )|| 2 + e in - (Vd p (A fc ), A)) . 

Taking now A = 0, k = fc ou t and using an inductive argument, we obtain: 

F(i i fcout ) - F* < Ld J A ° 112 + e in = 1, (33) 

fcout 4 

provided that A 0 = 0. From (l30l) . (l3ll) and m, we obtain that the average 
primal iterate u kout is e-primal optimal. □ 


Further, we analyze the primal convergence rate of algorithm DFGM in the 
average primal iterate u k : 

Theorem 4 Let e > 0 be some desired accuracy and u k be the primal av¬ 
erage iterate given in generated by algorithm DFGM, i.e. Algorithm 

DFOM(d p ,lCd) with 6k +1 = 1+ ^ +4 — for all k > 0 , using the inner accu¬ 
racy from USD- Then, after k ou t number of outer iterations given in m , u k e ° ut 
is e-primal optimal for the original QP CD- 


k 

Proof Recall that we have defined Sk = X) 6k- Then, it follows: 

j =o 


—< 6k < k and Sk = 6\. (34) 

For any j > 0 we denote z J = /i- 7 + ^j^Vdpipi) and thus we have A- 7 = [V] . 

In these settings, we have the following relations: 


(J-V 


\2L d 


Vd p (^) - {z> 



= 6 , 


^ + 2L^ dp ^ 


-A - 7 


>Cd 


= 0j(X j - ij?) 

= 0j{X j - A 7-1 ) + (0,-1 - 1)(A - 7-2 - A- 7-1 ) 

= A - 7-1 + 9j(\ j - A- 7-1 ) - (A - 7 " 2 + 6^-1 (A - 7 " 1 - 

'-—^ --- 


A' 7-2 )) • (35) 
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For simplicity consider A 2 = A 1 = A 0 and 9-i = 9 d = 0. Adding up the 
above equality for j = 0 to j = fc, multiplying by and observing that 
s(u J , /+) + H — [zi]ic d £ /C for all j > 0 , we obtain: 


dist*; ( Gu k + g) < 


it ( GuJ+g ~ s (V>h) - 2i d(~ j - [*»y 


J=0 

fc 






# _ i°\\ < , 4Ld l|l fc - H|. 


3=0 

Ld 

S k 


(fc + 1) 2 


Taking A = A* in Lemma [3] and using that the two terms 9^,(F* — d p (X k )) and 

k 

9jA(X*,fjP) are positive, we get: 

3=0 


I ^ 


A*||< a ||A 0 -A *|| 2 + ^TlA 1 < ||A°- 

i= 1 

< HA 0 - Al + (|g ) 1/2 (fc + 1)3/2 


*'■ + \lk {k+v 

Vfc > 0. 


Thus, we can further bound the primal infeasibility as follows: 

I,, k ,on „ 4Td 


dist^; ( Gu k + g) < 


-HI < 


< 


(fc + 1) 2 " "-(fc + 1) 2 

8 id-Rd , 0 f LdCjnX 1 ^ 2 


— A*|| + Rd) 


(fc + 1) 2 \k + lj ■ 

Therefore, using fc out and ei n from (I18fl . it can be derived that: 


distAc(Gu^ out +g) < 


SL d R d (L d ei n \ 1/2 3e 


^ut 


\ fcout / 


< 


Rd 


(36) 


(37) 


Further, we derive sublinear estimates for primal suboptimality. First, note 
the following relations: 

A(\,g k ) = d p (g k ) + (Vd p (/i fc ), A - /) - d p ( A) 

= C P (u k , g k ) + <Vd p (/), A - ft k ) - dp(X) 

= F(u k ) + (A, Gu k + g — s(u k ,g k )) + ||| Gu k +g- s(u k ,g k )|| 2 - d p ( A) 
> min F(u k ) + (A, Gu k + g — s) + ^||Gw fc + g — s|| 2 — d p { A) 

sG/C 2 

= £ p (u k ,X)-d p (X). 
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Summing on the history and using the convexity of £p(*, A), we get: 
k k 

> ^(£ p (zF, A) - d p ( A)) 

3=0 j=l 

> S k (. C p (u k e , A) - d p ( A)) = ^ A) - d p (A)) . (38) 

Using (1551) in LemmajSJ and dropping the term L d \\l k — A|| 2 , we have: 


F{u k ) + {Gu k + g — s{u k , A),A) 


dp(A fe )<^||A° 

a k 


Af + 


(39) 


Moreover, we have that: 



Now, by choosing the Lagrange multiplier A = 0 and k = fc ou t in (I551l . we have: 
F{u k e ° M ) - F* < F(u kout ) - d p ( A fcout ) < 2Ld 2 Ri + 2fcoutei„ < (40) 

iDllt 4 


On the other hand, we have: 

F* = min F(u) + (\*,Gu + g-s)<F{u k ) + (\*,Gu k +g~\Gu k +g] r ) 
< F{u k ) + FL d dist K (GMe + g). 

Taking k = fc 0 ut and £i n from (fT5|) . and using (1571) . we obtain: 

-3 e< F{u k ° M ) ~ F*. (41) 

Finally, from (1571) . (HOI) and (TIlT) . we get that the primal average sequence u kout 
is e primal optimal. □ 

In conclusion, in DuQuad we generate two approximate primal solutions u k 
and u k for each algorithm DGM and DFGM. From previous discussion it 
can be seen that theoretically, the average of primal iterates sequence u k has a 
better behavior than the last iterate sequence u k . On the other hand, from our 
practical experience (see also Section ED we have observed that usually dual 
first order methods are converging faster in the primal last iterate than in a 
primal average sequence. Moreover, from our unified analysis we can conclude 
that for both approaches, ordinary dual with Q 0 and augmented dual with 
Q y 0, the rates of convergence of algorithm DFOM are the same. 
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5 Total computational complexity in DuQuad 


In this section we derive the total computational complexity of the algorithmic 
framework DFOM. Without lose of generality, we make the assumptions: 
Rd > 1, e < 1 , A max (Q) > ||G|| 2 . However, if any of these assumptions does not 
hold, then our result are still valid with minor changes in constants. Now, we 
are ready to derive the total number of iterations for DFOM, i.e. the total 
number of projections on the set U and of matrix-vector multiplications Qu 
and G T A. 


Theorem 5 Let e > 0 be some desired accuracy and the inner accuracy ei n 

8 ff? 

and the number of outer iterations k ou t be as in (1181) . By setting p = 
and assuming that the primal iterate u k is obtained by running the Algorithm 
FOM(C p {• : p k ),U), then u k (u k ) is y/e (e) primal optimal after a total num¬ 
ber of projections on the set U and of matrix-vector multiplications Qu and 
G T A given by: 


ktotal 


2±\\G\\D u R d 

e 

16||G||J? d , f S\\G\\DuR d 
y/aZe ° \ e 


if = 0 

if > 0 . 


Proof From Lemma Q] we have that the inner problem (i.e. finding the pri¬ 
mal iterate u k ) for a given p k can be solved in sublinear (linear) time using 
Algorithm FOM(£ p (^ 1 •), U), provided that the inner problem has smooth 
(strongly) convex objective function, i.e. C p {p k , •) has ac = 0 {pc. > 0). More 
precisely, from Lemma [T] it follows that, regardless if we apply algorithms 
DFGM or DGM, we need to perform the following number of inner itera¬ 
tions for finding the primal iterate u k for a given pf: 



if (t c -— 0 
if <Jc > 0 . 


Combining these estimates with the expressions (1181) for the inner accuracy ei n , 
we obtain, in the first case ac = 0 , the following inner complexity estimates: 


I'h 


( 8 LzDj, 


1/2 


pL c D? r ) 1 ' 2 (2L d R 2 d ) 1 ' i 

e 3/4 i 


if DGM 
if DFGM. 


Multiplying k- ln with the number of outer iterations A: out from (1181) and mini¬ 
mizing the product k- m k on t over the smoothing parameter p (recall that Lc = 
Amax(Q) + p||G || 2 and L& = a - (QH-p||G || 2 we °bt a i n the following optimal 
computational complexity estimate (number of projections on the set U and 
evaluations of Qu and G T A): 


^total — {kout kin)* — 


2A\\G\\DuR d 


e 
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which is attained for the optimal parameter choice: 



e 


Using the same reasoning for the second case when ac > 0, we observe that 
8R 2 

the value p = — rl is also optimal for this case in the following sense: the 

difference between the estimates obtained with the exact optimal p and the 
8R 2 

value —- are only minor changes in constants. Therefore, when ac > 0, 
the total computational complexity (number of projections on the set U and 
evaluations of Qu and G T A) is: 


*t*otal = (fcoutfcin)* 


16||G||fld ( 8\\G\\DuRd 

\[vc£ 8 V e 


□ 

In conclusion, the last primal iterate u\ is ^/e-primal optimal after O(-) 
(<D(-^ logy)) total number of projections on the set U and of matrix-vector 

multiplications Qu and G T A, provided that ac = 0 (ac > 0). Similarly, the 
average of primal iterate u * is e-primal optimal after 0(\) (0(-^ log 7 )) total 
number of projections on the set U and of matrix-vector multiplications Qu 
and G T A, provided that ac = 0 (ac > 0). Moreover, the optimal choice for 
the parameter p is of order 0( 7 ), provided that A m i n (Q) = 0. 


5.1 What is the main computational bottleneck in DuQuad? 

Let us analyze now the computational cost per inner and outer iteration for 
Algorithm DFOM(d p ,/C*) for solving approximately the original QP ([Tjl: 

Inner iteration: When we solve the inner problem with the Nesterov’s algo¬ 
rithm FOM(£ p (/x, •), C7), the main computational effort is done in computing 
the gradient of the augmented Lagrangian C p (p,-) defined in ([3]), which e.g. 
has the form: 


V/2p(/x, u) — (Q + pG T G)u + (q + G T p + pG T g). 

In DuQuad these matrix-vector operations are implemented efficiently in C 
(matrices that do not change along iterations are computed once and only 
G T p is computed at each outer iteration). The cost for computing VC p (p 1 u) 
for general QPs is 0(n 2 ). However, when the matrices Q and G are sparse 
(e.g. network utility maximization problem) the cost 0(n 2 ) can be reduced 
substantially. The other operations in Algorithm FOM (C p (p,-),U) are just 
vector operations and thus they are of order 0(n). Thus, the dominant oper¬ 
ation at the inner stage is the matrix-vector product. 
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Outer iteration: When solving the outer (dual) problem with Algorithm 
DFOM(d p , 1C*), the main computational effort is done in computing the in¬ 
exact gradient of the dual function: 

Vt ip(ii) = Gu(p) + g- s(zi(/x), p). 

The cost for computing Vd p {p) for general QPs is 0{np). However, when the 
matrix G is sparse, this cost can be reduced. The other operations in Algorithm 
DFOM(d p , ICd) are of order 0{p). Thus the dominant operation at the outer 
stage is also the matrix-vector product. 

Fig. [2] displays the result of profiling the code with gprof. In this simulation, a 
standard QP with inequality constraints and dimensions n = 150 and p = 225 
was solved by Algorithm DFGM. The profiling summary is listed in the order 
of the time spent in each file. This figure shows that almost all the time 
for executing the program is spent in the library module math-functions, c. 
Furthermore, mtx-vec-mul is by far the dominating function in this list. This 
function is multiplying a matrix with a vector, which is defined as a special 
type of matrix multiplication. 


Name (location) 

Sampl Calls 

Time/Call 

%Time 

T Summary 

154 



100.0% 

► dfgm.c 

0 



0.0% 

► fgm.c 

1 



*65% 

► generaljunctions.c 

0 



0.0% 

► init_problem.c 

0 



0.0% 

► main.c 

0 



0.0% 

T mathjunctions.c 

153 



99.35% 

abs_2 

0 

12983 

0ns 

0.0% 

mtxjranspose 

0 

1 

0ns 

0.0% 


n 

92.21% 

obj 

0 

13234 

0ns 

0.0% 

► vector_add 

2 

26464 

755ns 

1.3* 

► vector_copy 

1 

27464 

364ns 

$.65% 

vectorjnax 

0 

12860 

0ns 

0.0% 

vector max with zero 

0 

372 

0ns 

0.0% 

► vector_min 

1 

12860 

777ns 

0.65% 

► vector_mul 

2 

26718 

748ns 

*3% 

vector_norm_2 

0 

124 

0ns 

0.0% 

► vector_scalar_mul 

3 

26215 

1.144us 

1.95% 

► vectorjub 

2 

26960 

741ns 

1.3% 


Fig. 2: Profiling the code with gprof. 


In conclusion, in DuQuad the main operations are the matrix-vector products. 
Therefore, DuQuad is adequate for solving QP problems on hardware with 
limited resources and capabilities, since it does not require any solver for linear 
systems or other complicating operations, while most of the existing solvers for 
QPs from the literature implementing e.g. active set or interior point methods 
require the capability of solving linear systems. On the other hand, DuQuad 
can be also used for solving large-scale sparse QP problems since the iterations 
are very cheap in this case (only sparse matrix-vector products). 
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6 Numerical simulations 

DuQuad is mainly intended for small to medium size, dense QP problems, but 
it is of course also possible to use DuQuad to solve (sparse) QP instances of 
large dimension. 


6.1 Distribution of DuQuad 

The DuQuad software package is available for download from: 

http: / / acse.pub.ro / person/ion-necoara 

and distributed under general public license to allow linking against propri¬ 
etary codes. Proceed to the menu point “Software” to obtain a zipped archive 
of the most current version of DuQuad. The users manual and extensive source 
code documentation are available here as well. 

An overview of the workflow in DuQuad is illustrated in Fig. [3] A QP problem 
is constructed using a Matlab script called test.m. Then, the function duquad.m 
is called with the problem as input and is regarded as a prepossessing stage 
for the online optimization. The binary MEX file is called, with the original 
problem and the extra info as input. The main.c file of the C-code includes 
the MEX framework and is able to convert the MATLAB data into C format. 
Furthermore, the converted data gets bundled into a C struct and passed as 
input to the algorithm that solves the problem. 


MaHab 


testm *K*i*ljn 



Fig. 3: DuQuad workflow. 
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6.2 Numerical tests: case K, = {0} 

We plot in Fig. 0] the average CPU time for several solvers, obtained by solv¬ 
ing 50 random QP’s with equality constraints (Q >: 0 and K, = {0}) for each 
dimension n, with an accuracy e = 0.01 and the stopping criteria | F(u*) — F 1 *! 
and 11 Gu 1 / +g \| less than the accuracy e. In both algorithms DGM and DFGM 
we consider the average of iterates u 1 /. Since Q >z 0, we have chosen p = 0(1/e). 
In the case of Algorithm DGM, at each outer iteration the inner problem is 
solved with accuracy ej n = e. For the Algorithm DFGM we consider two sce¬ 
narios: in the first one, the inner problem is solved with accuracy ej n = 0.001, 
while in the second one we use the theoretic inner accuracy m- We observe 
a good behavior of Algorithm DFGM, comparable to Cplex and Gurobi. 



Fig. 4: Average CPU time (ms) for solving QPs, with Q > 0 and K = {0} of 
different dimensions, with several solvers. 


6.3 Numerical tests: case K, = 

We plot in Fig. [5] the number of iterations of Algorithms DGM and DFGM 
in the primal last and average iterates for 25 random QPs with inequality 
constraints (Q >*- 0 and /C = M0) of variable dimension ranging from n = 10 
to n = 500. We choose the accuracy e = 0.01 and the stopping criteria was 
| F(u) — F*| and dist;t(GM + g ) less than the accuracy e. From this figure we 
observe that the number of iterations are not varying much for different test 
cases and also that the number of iterations are mildly dependent on problem’s 
dimension. Finally, we observe that dual first order methods perform usually 
better in the primal last iterate than in the average of primal iterates. 
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Fig. 5: Number of outer iterations on random QPs (Q >~ 0, KL = R0) for DGM 
and DFGM in primal last/average of iterates for different test cases of the 
same dimension (left) and of variable dimension (right). 
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